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Abstract 

Two fourth order accurate Osher schemes are pre- 
sented which maintain higher order accuracy on non- 
uniform grids. They use either a conservative finite dif- 
ference or finite volume discretization. Both methods 
axe successfully used for direct numerical simulations of 
flat plate boundary layer instability at different Mach 
numbers. Results of growth rates of Tollmien- Schlicht- 
ing waves compare well with direct simulations of in- 
compressible flow and for compressible flow with results 
obtained by solving the parabolic stability equations. 

1. Introduction 

The study of boundary layer stability and transition to 
turbulence is one of the classical topics in fluid mechan- 
ics. Linear and weakly non-linear theory, together with 
experiments, have been successful in describing several 
of the important instability mechanisms in compress- 
ible boundary layers, e.g. Tollmien- Schlichting (TS) 
waves, or first modes, and higher modes, which come 
into play at supersonic Mach numbers, Mack 15 , and 
theory was also successful in describing secondary in- 
stability, Herbert 10 . 

The complicated phenomena in the non-linear stages, 
leading to transition and turbulence, however, require 
further understanding. Direct numerical simulations 
can provide some of this information, but their appli- 
cation to compressible boundary layers has been hin- 
dered by many obstacles. To mention just a few, high 
order accuracy is required on non-uniform grids and a 
severe time step limitation is encountered due to the 
small grid spacing in the boundary layer when using an 
explicit time integration method; whereas with implicit 
time integration methods it is difficult to maintain time 
accuracy. 
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Spectral methods have been very successful in simu- 
lating incompressible flows in simple geometries, such 
as channel flow, e.g. Laurien and Kleiser 11 , but are 
not easily extended to more complicated geometries. 
The recently popular compact finite difference schemes, 
Lele 13 , do not have the geometric limitations of spectral 
methods and have been successfully applied to mixing 
layers and shock turbulence interaction, Lee et al 14 . 
Unfortunately, compact finite difference schemes and 
also spectral methods, cannot capture shocks and if 
they appear they have to be fully resolved, which can 
require prohibitively small grid spacings, Lee et al. 14 . 

There have been several attempts to use finite differ- 
ence schemes for direct simulation of transition in com- 
pressible boundary layers. The most frequently used 
method is the fourth order accurate version of Mac Cor- 
mack’s scheme, developed by Gottlieb and Turkel 7 , e.g. 
Maestrello et al. 16 and Bestek et al. 3 . This method can 
only achieve higher order accuracy on grids generated 
as the product of two or more one-dimensional analytic 
transformations, limiting its applicability to relatively 
simple, smooth flows. 

Most frequently explicit time integration methods have 
been used, but for many transitional flows the Courant- 
Friedrichs-Lewy (CFL) time step limitation is not nec- 
essary to maintain time accuracy. Recently Rai and 
Moin 21 developed a numerical scheme which solves the 
compressible Navier-Stokes equations using a time ac- 
curate upwind biased implicit method and were able 
to simulate bypass transition. This method alleviates 
the time step limitation of explicit methods, but has 
as main drawback that it uses the non-conservative 
form of the Navier-Stokes equations and only allows 
grid stretching in one direction. The grid stretching, 
however, does not have to be analytically defined be- 
cause the higher order finite difference approximations 
are generated numerically in physical space. 

In this paper two alternative methods will be dis- 
cussed. The first method is a higher order accurate 
extension of the MUSCL scheme, originally developed 
by Van Leer 27 as a second order accurate extension of 
Godunov schemes. The scheme is related to the multi- 
dimensional essentially non-oscillatory (ENO) schemes 
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developed by Casper and Atkins 4 and Harten et al. 9 . 
The second method is a higher order accurate upwind 
biased version of the Osher scheme, which main tains 
its high order accuracy on non-uniform grids. Higher 
order accurate Osher schemes were also discussed by 
Rai 20 , but his method is only higher order accurate on 
a uniform grid. 

The discussion in this paper will be restricted to smooth 
flows, but both schemes have been extensively tested 
for flows with shocks and other discontinuities in one- 
dimensional flows. A detailed discussion of the benefits 
and application of these schemes to non-smooth flows 
can be found in Van Der Vegt 25 . The purpose of this 
paper is to investigate if these methods are accurate 
and efficient enough to be used as tools for direct sim- 
ulation of boundary layer instability and transition to 
turbulence. 

The finite volume scheme has as main benefits that it 
is a truly multi-dimensional scheme, whereas the finite 
difference scheme uses dimensional splitting. The finite 
volume scheme automatically satisfies conservation and 
is most closely related to the integral formulation of 
the compressible Euler and Navier-Stokes equations. It 
also maintains higher order accuracy at sonic points, 
which is not true for the finite difference formulation. 
The finite volume method, however, is slightly more 
costly than the finite difference scheme and requires 
significantly more effort to implement. 

The use of an upwind scheme is beneficial for direct 
simulations of compressible flow, because it automat- 
ically controls aliasing errors. There are two types 
of upwind schemes, those based on flux vector split- 
ting, e.g. Van Leer 28 and St eger- Warming 24 , and those 
based on a Godunov approach which use the solution 
of a Riemann problem. Godunov schemes most closely 
mimic the physics of wave propagation in compressible 
flow and have excellent shock capturing properties. In 
this class of schemes the Roe and Osher approximate 
Riemann solvers are the most popular, see Roe 22 and 
Osher and Solomon 18 . The Osher scheme has been cho- 
sen because it has a very low numerical dissipation in 
boundary layers, which is crucial for direct simulations, 
and it has a continuously differentiable flux, which is 
important for implicit schemes. In addition it satisfies 
the entropy condition and has good shock capturing 
properties. 

In the next two sections of this paper first the basic 
equations will be discussed and the higher order accu- 
rate numerical schemes will be presented. The paper 
will conclude with a discussion of the results of compu- 
tations of boundary layer instability at various Mach 
numbers. 


2. Navier-Stokes Equations 

The Navier-Stokes equations can be considered either 
in integral formulation, leading to the finite volume dis- 
cretization, or in differential form, which is the basis for 
the finite difference discretization. The finite volume 
method automatically satisfies the conservation prop- 
erties of the equations but care has to be taken that the 
finite difference method is in the so-called conservation 
form. It is otherwise not possible to obtain a weak so- 
lution with the proper jump relations at discontinuities 
in the limiting case of vanishing viscosity, as demon- 
strated by Lax and Wendroff 12 . A detailed discussion 
of finite volume and finite difference methods and their 
differences can be found in Vinokur 29 . 

The integral formulation of the compressible Navier- 
Stokes equations is defined as: 

[ U dV- [ U dV+ [ i T(V)ndSdt = 0 
Jn(t 2 ) J n(t!) Ju Jdn{t) 

(2- 1 ) 

Here Q(f) is the flow domain with boundary 90(f) at 
time f and n the unit outward normal vector at 90. 
The vector U represents the conserved variables: 
(p y pu,pv y pw y e ) T , with p density, u = (u y v y w) T the 
Cartesian velocity components, and e total energy. The 
matrix T, which represents the fluxes through the sur- 
face 90(f), consists of two parts, F the inviscid flux, 
and V the viscous flux, with T—Y-V . The inviscid flux 
contribution F has as components: 


Fi 
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p uv 
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where the dimensionless pressure p is determined from 
the equation of state: p = pc v T/(~fM 2 ), with c*, the 
specific heat at constant volume, 7 the ratio of specific 
heat at constant pressure and constant volume, M the 
Mach number and T temperature. The temperature T 
is given by the relation: T = 7(7 — l)M 2 (e — |pu ■ 
u)/{pc v ). 


2 



The viscous contribution V has as components: 
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with the stress tensor r and the variables /3 X , j3 y and 
f3 z defined as: 


( . *,dv dw \ . 

r xx = [i2fi + \)- + \(-+—)j/Re 

( ,du dv 
= /Re 

( ,du dw.\ . 

^=( (2 ' ,+ < +A( £ + £ ) ) /j& 

*=(*S + b>)/- 

k dT 

( 3 X =ur xx + vr xy + wt X2 + (7 _ 1)M - 2p - ^ 

K dT e . 

^ =UT Xy + VTyy + IDT,, + (t _ ^ ( 2 ‘ 8 ) 

is dT 

A =.r„ + w„ + «r„ + (t _ 1)AfIf , r Tz 


Here Re represents the Reynolds number, Pr Prandtl 
number, k the coefficient of thermal conductivity and 
fi and A the first and second viscosity coefficient. All 
computations were done using the relation A = — |m? 
with & given by Sutherland’s law. The non-dimensional 
variables are defined with respect to the freest ream ve- 
locity, density, temperature, viscosity, thermal conduc- 
tivity and specific heat. 


If we assume that all variables are continuously differ- 
entiable in time it is possible to rewrite equation (2.1) 
into: 


f U dV + <f T(U) ■ ndS = 0 (2.6) 

at Jn(t) Jdn{t ) 

A special constraint can be derived from this expres- 
sion, namely the geometric conservation law. Inserting 
a uniform flow field in equation (2.6) we obtain: 

/ n dS = 0 (2-7) 

Jan(t) 


which states that the cell face dH(f) must be closed. 
When dividing the total flow field in a set of non- 
overlapping cells this constraint puts limitations on the 
way how to compute the cell faces and volumes. They 
all have to add up to the total volume and each cell 
will have to be closed otherwise a uniform flow field 
will be disturbed due to contributions from the metrics. 
This is a non-trivial problem when deriving higher or- 
der schemes and will be discussed in the next sections. 

The differential form of the compressible Navier-Stokes 
equations is directly obtained from the integral formu- 
lation, equation (2.6), using Gauss’ theorem and con- 
sidering an arbitrary volume 0: 

^+V-T(U) = 0 (2.8) 

This is the conservative formulation of the compress- 
ible Navier-Stokes equations and it is important that 
the discretization also can be written in conservation 
form. The geometric conservation law must be satis- 
fied and this requires great care in dealing with the 
metrical coefficients in a finite difference discretization, 
especially in three dimensions. 

3. Higher Order Accurate Finite Volume 
Scheme 

Second order accurate finite volume schemes have been 
around for a long time. When extending the accuracy 
beyond second order several problems are encountered. 
It is no longer valid to equate cell averaged values with 
values at the cell center and a more elaborate algo- 
rithm to reconstruct point values from cell averaged 
values is required. In addition one has to compute the 
integrals of the fluxes over the cell surfaces more ac- 
curately. The standard procedure of multiplying the 
flux with the cell surface is only second order accurate. 
The geometry in a higher order accurate finite volume 
method also has to be represented more accurately, es- 
pecially at the boundary. The relations for cell area 
and volume as given by Vinokur 29 are no longer suf- 
ficient. They are exact for hexahedrons with straight 
line edges, but not for cells with curved edges. 

Until now very few attempts have been made to de- 
velop higher order accurate finite volume methods. For 
structured grids Qarten et al. 9 and Casper and Atkins 4 
developed the multi-dimensional ENO schemes and 
Abgrall 1 and Harten and Chakravarthy* studied these 
schemes for unstructered grids. Despite the fact that 
the unstructured approach has more flexibility in treat- 
ing complicated geometries and allows local grid re- 
finement it was decided to use the structured grid ap- 
proach developed by Casper and Atkins 4 . Both CPU 
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time and memory usage in the unstiuctered schemes, 
is so large that it is not feasible to use these meth- 
ods for direct or large eddy simulations of compressible 
flow. The structured grid approach, however, requires 
a smooth C 3 grid. It is possible to deal with locally 
non-smooth boundaries and intersections in the finite 
volume approach by subdividing the grid into smooth 
subsections using a multi- block approach. The smooth- 
ness requirements of the grid will limit the application 
of this method to fairly simple geometries, but due to 
their high cost large eddy and direct simulations will 
be limited to simple flows for quite a while. The finite 
volume method is, however, considerably more flexible 
than spectral methods which require a C°° grid. 

The use of a structured grid makes it possible to trans- 
form the physical domain to a simpler computational 
domain. Let £, 77 , and £ be the coordinates in the com- 
putational domain then they are related to the physical 
coordinates x y y, z by the relations: 

i =£(*, 3 l,z) 

i]=7){x,y,z) (3.1) 

< =<(z,y,z) 

The finite volume discretization on a structured grid 
is obtained by dividing the domain ft into a regular 
partition ftt f j f *. Each element is a hexahedron 

with coordinates s*, yj and z\ \ for the top right corner. 
The subdomains are non-overlapping and their 

sum is equal to the domain ft. 

The cell average in a cell with index (i, j, k) is 

defined as: 


U. 


*»j»* — 


Fo/(ft. 


i [ 


VdV 


(3.2) 


with Volfoij^k) the volume of cell fttj,*. The equation 
for the cell average is obtained by limiting the integra- 
tion domain ft to ft*,j,k in equation (2.6), and is equal 
to: 

^U i , i , Jb (<)+h i , J , t (U) = 0 (3.3) 

with the flux integral hi t j t k at the surface dftij^ of the 
cell fttj.jk defined as: 


hi.i.JbW = 


1 [ 

i f 


JF(U) ■ ndS 


F(U )dS 


(3.4) 


” Voi(ft* 

with F the flux normal to the cell surface: 
F = kjFi + k 2 F 2 4* 


(3.5) 


and n = (ki, k 2 > ks) T . This relation gives the method of 
lines formulation for the cell averaged equation Ui lJ} *(i). 
The numerical flux in a higher order finite volume scheme 
now is constructed such that it approximates the exact 
flux at time t = [n + l)At up to 0{h T ): 

USI - [Eh( At) • Tf] . _ k = 0(h r ) (3.6) 

with Eh(&t) the numerical solution operator. The 
higher order accurate finite volume scheme therefore 
gives an r-th order accurate approximation to the cell 
averages, not the point values as in a finite difference 
scheme. 

The two important ingredients of a higher order accu- 
rate finite volume method are the reconstruction of the 
point values U(x) from the cell averaged values Ui,j,k 
and the computation of the flux F at the cell face. The 
point values are necessary to compute the fluxes 
at the cell faces which depend on U(x) and not on 
U. This is done with the reconstruction method dis- 
cussed in section 3.1. The fluxes at the cell faces are 
computed by considering a Riemann problem at each 
cell face. This method was introduced by Godunov 6 
for first order accurate schemes and extended by Van 
Leer 27 to second order accuracy. Instead of using the 
exact solution of the Riemann problem the approximate 
Riemann solver developed by Osher and Solomon 18 is 
used, because it is less expensive than the exact solu- 
tion, but has excellent shock capturing properties and 
minimal dissipation in boundary layers. 

3.1 Reconstruction of Point Values from, Cell Averages 

In one dimension the most successful reconstruction 
technique is based on the primitive function method, 
see e.g. Harten et al. 9 . This method was extended by 
Harten et al. 9 and Casper and Atkins 4 to multiple di- 
mensions using the one-dimensional primitive function 
reconstruction technique first to compute line averages 
from cell averages, after which point values are com- 
puted with a second one-dimensional reconstruction. It 
is, however, possible to define a primitive function di- 
rectly for the multi-dimensional problem without hav- 
ing to use one-dimensional primitive functions. 

Define the primitive function U as: 

W(*,1 J,0= f f f* 

" £o drjo * £0 

(3.7) 

then the flow field U satisfies the relation 

p - 8 > 
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with £, 7j and £ the coordinates in computational space 
and J the Jacobian of the coordinate transformation. 
The primitive function U can be related to the cell av- 
erage U in physical space using the following relation: 

M{£iiVjyCk) = 


This relation is the basis for the higher order finite vol- 
ume scheme. It defines the primitive function U di- 
rectly at the corners of each hexahedron in computa- 
tional space (CiyVji 00 and is conservative. The point 
values U(x) are then obtained using equation (3.8). 

When the flow field is smooth the following procedure 
can be used to compute the pointwise data: For a sur- 
face in the plane £ = £*, first approximate the left and 
right states by differentiating the primitive function U 
with upwind biased differences in both f -directions with 
fourth order accuracy for all indices j and k. Then the r\ 
and £-derivates are computed at the Gauss quadrature 
points and divided by the local Jacobian to obtain the 
point values. This process is repeated for the planes 
with Tf — tj* and C — C* and works well for smooth 
flows. For flows with discontinuities the ENO recon- 
struction, which tries to use only data from the smooth 
part of the flow field has to be used. This will be dis- 
cussed in a future paper. 

Although this process is rather simple care has to be 
taken to prevent loss of accuracy due to truncation er- 
rors, because the primitive function U frequently be- 
comes very large or small. The way to prevent this is to 
construct the primitive function using only those cells 
around the cell with index (i, j, A ) which are needed to 
compute the derivatives in equation (3.8) and dynami- 
cally adapt the indices i 0 > Jo and k$. Accuracy is further 
improved by combining the process of summation and 
differentiation, e.g. first compute the sum with i-index, 
then differentiate this result and compute the summa- 
tion with j-index, and so on. 

In order to preserve uniform flow it is necessary to 
compute the Jacobian of the coordinate transformation 
at the Gauss quadrature points the same way as done 
for the flow field U. This can be accomplished most 


/*u rvj rii 

/ / / U\J\d$'dT)'dC 

•'Co Jfo 

k j i k> jt 

SEE/ j \J\J\di'dv'dC 

E E E / v(*,t)dv 

k'=k 0 j’=joi'=io 
k j i __ 

E E E v<*(niw)*iw 

k'=:k o j'=jo t f =»0 




easily by multiplying equation (3.8) with the Jacobian 
and inserting a uniform flow field in equations (3.8-9) 
will give the Jacobian at the Gauss quadrature points 
with the same reconstruction process as for U. This 
procedure is important because otherwise the recon- 
struction process will generate small errors which can 
be very annoying on stretched grids. 

3.2 Higher Order Accurate Flux Integrals 

The discretization of the integral formulation of the 
compressible Navier-Stokes equations (3.3) is completed 
with the approximation of the flux integrals given by 
equation (3.4). Casper and Atkins 4 gave a detailed 
analysis of the accuracy required in the reconstruction 
problem and the computation of the flux integrals to 
obtain an accurate solution with an error of order h T . 
In this paper we limit the discussion to fourth order ac- 
curacy. The use of a Gauss quadrature method is the 
most efficient way to compute the flux integrals with 
fourth order accuracy. First the integral over the cell 
boundary is split into the sum of integrals over the six 
cell faces: 


Km(u) = 




i.L 


F(U)dS (3.10) 


here the index l refers to one of the six faces of the 
hexahedron. Then the flux integrals at each cell face 
can be further evaluated using the Gauss quadrature 
rule at each cell face and the Osher approximate Rie- 
mann solver 17 * 18 is used to compute the flux at each 
quadrature point: 


J F(U)dS = 




=\ J (f l +F r - j \8F\dU j \J\dS 
aft!,. ^ r / 


F L + F R - J\dF\dU\ \J 3m ,,\dS 

r ' ft.,1 

(3.11) 

with , the cell face with index l in computational 
space. The indices g^i refer to the quadrature points 
with index m in the cell face where the fluxes of the 
left and right states and the Osher path integral 

are computed. The quadrature points have coordinates 

2 ± awiuniag that the hexahedrons 

sides have unit length. This relation is used by Casper 
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and Atkins 4 and requires four calculations of the ap- 
proximate Riemann flux for each cell face which sig- 
nificantly increases the computing time, Harten and 
Chakravarthy 8 suggested that only one computation of 
the (approximate) Riemann flux is required to main- 
tain accuracy in smooth flows. Due to the fact that the 
Riemann flux, and also the approximate Riemann flux 
according to Osher, is Lipschitz continuous and |Ux, — 
U R \ = O(JT) in smooth flows it is easy to show that the 
third component in the integral, equation (3. 11), can be 
approximated as: 

/ [ \9V\dU\J\dS = f |3F Ct \dU\J cl \ + 0(h r ) 

Jan'. h Jr Jr 

(3- 12 ) 

Here the index c/ refers to center of the cell face with 
index l. This relation significantly reduces the com- 
puting time, while maintaining the total accuracy at 
the slight expense of computing Ul.h at the cell face 
center* In regions with discontinuities it is, however, 
advisable to use the (approximate) Riemann solution 
at all the Gauss quadrature points. 

4. Higher Order Accurate Finite Difference 

Scheme 

The most difficult problem in deriving a higher or- 
der accurate finite difference scheme is to And a way 
to maintain that accuracy on a non-uniform grid. In 
upwind finite difference schemes, either based on flux 
vector splitting or using approximate Riemann solvers, 
the flux is a function of more than one grid point. When 
deriving the expression for the higher order differences 
care has to be taken to account for the changing met- 
rics, but this is frequently neglected. For instance the 
higher order Osher scheme derived by Rai 20 is only 
higher order accurate on a uniform grid. In addition 
to the accuracy requirement care has to be taken that 
the scheme is in conservation form and maintains uni- 
form flow, which is a non-trivial requirement for a finite 
difference scheme. The conservation property is impor- 
tant on physical grounds, the equations express con- 
servation of mass, momentum and energy, but is also 
important when dealing with discontinuities. In this 
paper, however, only smooth flow fields will be consid- 
ered. 

Before deriving the higher order Osher scheme it is 
necessary to study the first order Osher scheme in more 
detail . The conservative approximation to using 
Osher’s scheme is defined sis: 


where the conservative flux is defined as: 

E;+l ^(Et+1 (£t+|) + 

f 3uE + (^ + i)dU+ [ aoB-te + j)rfU) 

JTi JTi 

(4 * 2) 

with equivalent relations for d v F and G. The symbol 
d\j represents partial differentiation with respect to U 
and for ease of notation the j and k indices are omitted. 
In this relation Si + 1 refers to the metrical coefficients 
which are computed at the point with index i+i.The 
integrals in equation (4.2) are computed along a path 
in phase space, I\, and using the fact that the Riemann 
invariants are constant along this path Osher was able 
to derive exact analytic expressions for these integrals, 
see 20,22 . It is important to note that although the path 
integral I\ is from i to t + 1 only metrical coefficients at 
one point must be used for consistency. As can be seen 
directly from equation (4.2), the flux E depends on 
the positions with indices t, i + | and i 4-1. Using 
a Taylor series expansion with respect to both i and 
i + 1 Rai 20 was able to derive higher order conserva- 
tive expressions for c^E*. The dependency on i -j- 
however, was neglected, which reduces the order of ac- 
curacy of the scheme on a non-uniform grid, even on 
mildly stretched grids The analytical derivations nec- 
essary to obtain Rai’s higher order Osher scheme are al- 
ready tedious and taking care of the changing metrical 
coefficients, which would require Taylor series expan- 
sions up to at least fourth order in three independent 
variables, becomes unwieldy. 

An alternative is to compute the higher order accurate 
flux approximations numerically. This is done using 
the flux-ENO scheme discussed in Van Der Vegt 25 , but 
the stencil switching, which is part of ENO schemes, is 
eliminated in this paper. As starting point a different 
formulation of the first order Osher scheme is used: 

dtE i= [ duE + (&-i)<ttJ+ [ 3uE-(&+i)dU 

J 3 JTi 

( 4 - 3 ) 

It is easy to show that both formulations are equivalent, 
see Osher and Solomon 18 . The higher order scheme is 
derived using a Newton interpolation to the fluxes. The 
Newton interpolation, however, cannot be directly ap- 
plied to the integrals in equation (4.3) because of the 
path integrals. In order to use the Newton interpola- 
tion we use some simple relations, which were derived 
by Shu and Osher 23 : The function f(x) can always ex- 
pressed as: 

f(x) = iiL J h{x ' )dx> 

* x r 




(4.1) 
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and using its primitive function: F(x) = J^ o h(x t )dx f 
the following relation is obtained: 

f(* i ) = P(ps i + ^)-F(* i -^) (4.4) 

These relations can be used to link the primitive func- 
tion F to the flux integrals dff : 

df + =f; + _ Ci 

=F+ [x i+ i , Xi _ j] - F + [* 4 _ • , *<_|] (4.5) 

=2F + [*i + i i x i-i > x >-|] 

dfr=C + i-C 

=F"[x <+ |,® i+ i] - (4.6) 

=2F-[* i+ |,* i+ j,* i _|] 

with: 

df- = / duB-(^ +i )dU 

■/r 4 ’ 

dff = f 8 V E + fe_i)dU 

and F[®i + *, . . . , x,] the k-th divided difference defined 
as: 


The higher order approximation to d £ E* is obtained us- 
ing equation (4.4) and adding the positive and negative 
contributions in equations (4.7-8): 

a { Ei = d £ F(* i+ i)-3 { F(*i_i) (4.9) 

This relation is conservative, fourth order accurate, and 
maintains its higher order accuracy on a non-uniform 
grid because the change in metrical coefficients is prop- 
erly taken care off by means of the Newton interpola- 
tion. It satisfies the geometrical conservation law be- 
cause for each of the integrals df f, i E {i — 2,i + 2}, 
appearing in the divided differences, the metrics are 
chosen at indices * + i E {i — 2, t + 2}. The geometric 
conservation law then is automatically satisfied because 
for uniform flow each of the integrals gives a zero con- 
tribution independent of the metrical coefficients. The 
additional cost of computing the divided differences is 
negligible compared to the computation of the integrals 
df^ and the scheme is as efficient as a scheme with an- 
alytically derived coefficients. One additional remark 
must be made about the first divided difference in equa- 
tions (4.7-8). Their value is unknown, but not needed, 
because in equation (4.9) only their difference is used, 
which is exactly the first order contribution and given 
by equations (4.5-6). 


Ffct+jfcj • * • } — — (F[x*+Jk, . . . , X*+i] F[*i+Jfc-li * * *> *»]) 

In the derivation of the higher order accurate scheme 
the specific functional form of the functions f ^ and F* 
is not needed, only their divided differences. 

The primitive function F is now approximated with a 
fifth order Newton polynomial using the divided differ- 
ences defined in eq. (4.5-6). The higher order divided 
differences can be easily obtained by further extending 
the divided difference tables given by equations (4.5- 
6). The nodes in the Newton interpolation are chosen 
such that an upwind biased scheme is obtained. This 
relation is then differentiated at: 


3 f F+ (*i+i) =F + [x,+i,a: l _i] 

+ F+[x i+ 

+ 2F + [®j + | , • 


— 2F + [x i+ |, • 

' ' > x «-|] 

d £ F-(® i+ i) =F-[x i+ |,x j+ i] 

|-F-[« <+ 

+ 2F _ [x i+ i,- 


+ 2F~[x i+ i,- 



5. Implementation of Higher Order Schemes 

In order to obtain the high accuracy necessary for direct 
sin&ulations a large number of grid points is required. 
In order to efficiently run the program with such large 
grids a general three-dimensional multi-block code has 
been written. This gives more flexibility in managing 
the large memory requirements and it is easier to gen- 
erate grids for more complicated geometries. The in- 
viscid contribution in the compressible Navier-Stokes 
equations is discretized using the procedure described 
in the previous sections. The viscous terms are im- 
plemented using a fourth order accurate, conservative 
central differences for the finite difference scheme, but 
the viscous contribution is only second order accurate 
for the finite volume scheme. The development of a 
fourth order accurate viscous discretization for the fi- 
nite volume scheme is currently in progress. Second 
order accurate implicit time integration is used and in 
order to maintain time accuracy a Newton procedure is 
chosen, equivalent to the one used by Rai and Moin 21 
and Rai 30 . In this method the equations are discretized 
with all the fluxes computed at the new time level n+1. 
The resulting non-linear equations are linearized us- 
ing a Newton procedure and the large system of lin- 
ear equations is solved iteratively. In 20,21 the diagonal 
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form of approximate factorization according to Pulliam 
and Chaussee 19 was chosen, which only gives a crude 
approximation to the solution of the linear system. In 
this paper a more accurate iterative scheme based on 
the zebra line Gauss- Seidel method is used and for each 
step in the Newton procedure the system is solved up to 
machine accuracy. This method is used plane by plane 
and fits well into the Newton procedure. 

The Newton procedure makes it possible to have an 
implicit spatial discretization which is of lower accuracy 
than the explicit part, but because each time step the 
Newton procedure is iterated time accuracy is main- 
tained and the higher order spatial accuracy is pre- 
served. Usually four Newton iterations are sufficient, 
but for convergence of the Newton scheme it is nec- 
essary to solve the linear system of the implicit time 
integration with high accuracy. 

The Newton procedure requires the computation of 
the Jacobian of the flux vector, which is quite cumber- 
some to derive, especially for the viscous terms. Cur- 
rently, both the exact Jacobian of the Osher fluxes and 
the approximation using the Jacobian of the Steger- 
Warming 24 flux vector splitting are used in the im- 
plicit time integration. The exact Jacobian of the Osher 
fluxes improves the convergence rate compared with the 
approximate Jacobian. It is, however, approximately 
three times more expensive to compute the exact Ja- 
cobian and for most cases the computing time for both 
schemes is about equal. For steady flows, where the Ja- 
cobian has to be updated only after a certain number 
of time steps the exact Jacobian is more efficient. 

The boundary conditions for the Osher scheme are 
implemented using the procedure proposed by Osher 
and Chakravarty 17 . In this method a Riemann initial- 
boundary value problem is solved instead of an initial 
value problem, which is used in the interior of the do- 
main. This method is very robust and elegant, but a 
significant effort is required to derive all relations to 
compute the boundary fluxes and Jacobians. 

6. Results and Discussion 

In this section results will be presented of two simu- 
lations of ribbon induced boundary layer instability to 
demonstrate the ability of the two numerical schemes 
discussed in this paper to accurately predict boundary 
layer instability and their possible use for simulations of 
turbulent and transitional boundary layers. Although 
the results in this section are all two dimensional the 
full three-dimensional discretization was used. 

As a first step it is crucial to have extremely accurate 
solutions of the mean flow. In many previous stud- 
ies an analytically defined mean flow was used, but 


this becomes difficult for flows in more general geome- 
tries. The use of an analytically defined mean flow also 
has as disadvantage that it does not exactly satisfy the 
discretized equations and will generate numerical tran- 
sients. 

In Figures 1 and 2 the mean flow profiles obtained 
with the finite difference scheme discussed in section 
4 are plotted at 10 equally spaced stations with Re xy 
the Reynolds number based on the distance from the 
nose of the plate, ranging from 50.000 to 320.000 ver- 
sus the similarity parameter rf = *\/i2e, with Re the 
Reynolds number based on plate length, x the stream- 
wise distance from the nose of the plate and y the nor- 
mal distance. The freestream Mach number is .08. The 
dimensionless normal velocity in Figure 2 is defined as: 
v = Vy/Re x . At the inlet a boundary layer solution was 
specified and the grid was chosen to approximately fol- 
low the streamlines and generated with the GridGen2d 
package. The grid is non-orthogonal in the interior and 
the normal grid spacing increases downstream. In the 
same plot also the compressible Blasius solution is plot- 
ted and it is clear that all lines completely collapse for 
the mean flow stream wise component, Figure 1. The 
same is true for the normal velocity, Figure 2, except for 
the asymptotic value outside the boundary layer, which 
is slightly higher than the Blasius solution. This slight 
difference is correct because in the boundary layer solu- 
tion the small displacement of the boundary layer is not 
accounted for. The comparison between the theoreti- 
cal and numerical solutions is remarkable considering 
the high Reynolds number of the base flow. Especially 
the accurate solution of normal velocity component is 
noteworthy, because its value is much smaller than the 
stream wise component and more difficult to compute. 
Most tests of numerical schemes on a flat plate bound- 
ary layer use a Reynolds number which is considerably 
lower and only show the streamwise velocity compo- 
nent. Accurate boundary layer profiles were already 
obtained with 35 points in the normal direction, but 
the total grid consisted of 336 x 80 points to provide 
sufficient accuracy for the direct simulations, discussed 
in the next part. This demonstrates that the Osher 
scheme is considerably more accurate in boundary lay- 
ers than schemes based on flux vector splitting, Van 
Der Vegt 26 . 

The mean flow profiles obtained with the finite vol- 
ume scheme discussed in section 3 are similar, but the 
finite volume scheme turns out to be more sensitive to 
the smoothness of the grid on highly stretched meshes. 
Care has to be taken that the grid for the finite vol- 
ume scheme is at least three times differentiable. The 
sensitivity to the grid smoothness of the finite volume 
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scheme is caused by the fact that in the reconstruction 
process the cell averaged flow fleld in equation 

(3.9) is multiplied with the cell volume. The cell vol- 
ume changes much more rapidly than the grid spac- 
ing. The finite difference scheme is less sensitive to 
the grid, because it uses dimensional splitting and de- 
pends therefore only on the smoothness of the grid in 
each coordinate direction. The sensitivity to the grid 
smoothness of the finite volume scheme certainly needs 
further attention before this scheme can be used for 
more general applications. 

Another problem when using a high order accurate 
scheme to obtain steady solutions for high Reynolds 
number flows is the slow convergence to steady state. 
Due to the minimal amount of numerical dissipation it 
takes a significant amount of computing time to elim- 
inate all transients. Convergence to steady state was 
improved using an implicit scheme and CFL numbers 
between 1000 and 5000 were used to obtain steady re- 
sults with the fourth order accurate schemes. In or- 
der to further speed up convergence the computations 
were started with a first order accurate scheme till the 
residue was significantly reduced, followed by the fourth 
order scheme till machine accuracy was obtained. 

The first simulation of boundary layer instability is a 
comparison with the direct simulations of incompress- 
ible flow about a flat plate done by Fasel et al. 5 . All 
parameters were chosen as close as possible to the one 
used in their computations. The free stream Mach 
number was .08 and the Reynolds number based on 
flat plate length Re was 100.000/m. The plate started 
at x = .5 and ended at x = 3.2 and was extended 
with a buffer region with slowly increased grid spacing 
to x = 8. for the finite difference calculations and to 
x = 5. for the finite volume calculations. The plate has 
300 grid points in streamwise direction and the buffer 
region consisted of 36 points. The purpose of the buffer 
region was to damp out transients and thereby min- 
imizing reflections. Accurate non-reflecting boundary 
conditions for the compressible Navier- Stokes equations 
are not yet available. 

First a steady boundary layer solution was computed 
and the maximum pointwise value of the residual was 
less than 10“ 8 . The flow was then disturbed in a small 
region by periodic suction and blowing. The suction 
strip started at xi — .908 and ended at x 2 = 1.13. The 
amplitude is given by the relation: 

f v = Asin(x)(l - cos(x))sin(/3f) 


x 2 — 


The amplitude A for the computations is .0001. The 
parameter /3 is equal to 10, which results in a frequency 
parameter F = 100. Here the frequency parameter is 
defined as: F = /? x 10 6 /Re. 

The blowing and suction starts in the region were the 
boundary layer is linearly stable. This has as ben- 
efit that transients, which occur due to the suction 
and blowing, will damp out and a cleaner Tollmien- 
Schlichting wave is obtained. The time trace of all the 
flow variables along a line which corresponds to the po- 
sition of m axim um amplification was written to a file 
and Fourier analyzed. The Fourier analyzed signal was 
then used to compute the growth rate of the TS wave. 
Figure 3 shows the comparison of the growth rate —a* 
of the streamwise disturbances with the results of Fasel 
et al. 5 . A negative value of the growth rate means that 
the disturbance is growing. The large initial distur- 
bances are caused by the suction and blowing, but af- 
ter x = 1.4 the result from the finite difference scheme 
compares well with that from Fasel et al. 5 , which also 
agree with the theoretical non-parallel results of Gast- 
ner. The finite volume results are slightly less accurate 
than the finite difference results. This is partly due 
to the fact that the viscous contribution in the finite 
volume scheme is only second order accurate. The sec- 
ond simulation, a boundary layer at M = 0.5, which 
is at a considerably higher Reynolds number gives vir- 
tually identical results for both methods. The growth 
rates of the disturbances of the normal velocity com- 
ponent, Figure 4, also compare well with the results 
from Fasel et al. 5 , but the curves are shifted slightly 
upstream. This component becomes earlier unstable 
than in the simulation of Fasel et al. 5 . The growth rate 
of the normal velocity disturbances, however, strongly 
depends on the vertical position and further research is 
required to obtain more accurate results for this compo- 
nent. Contrary to the streamwise component there are 
no theoretical results for the growth rate of the normal 
velocity disturbances. 

The same procedure as used for the M = .08 boundary 
layer was used for a simulation of a flat plate boundary 
layer at Mach number M — .5. The parameters were 
chosen equal to the calculations done by Bertolotti 2 
using the linear Parabolic Stability Equations (PSE). 
This method takes the non-parallel effects of the bound- 
ary layer into account contrary to linear stability the- 
ory. The Reynolds number based on plate length was 
500.000/m. The simulations with the finite difference 
sch em e were done on two grids. The coarse grid has 
336 x 80 grid points and the plate started at x — .5 and 
ended at x = 3.2. Suction and blowing was started at 
x = .725 and ended at x ;= .86. The second grid con- 
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sisted of 436 x 80 points and the plate started at x = .5 
and ended at x — 2.3. Suction and blowing started at 
x = .5225 and ended at x = .6125. The buffer region 
ends at x = 8. for the coarse grid and at * = 5.1 for the 
fine grid. The first grid has about 10 grid points per 
TS wave and the fine grid has 20 point per wavelength. 
The frequency parameter F was equal to 20 and the 
free stream temperature Too — 20 6K. 

This case is more difficult than the incompressible sim- 
ulation because of the higher Reynolds number and 
many more TS wave periods have to be covered. The 
initial amplitude of the disturbances A was .0001 and 
Figure 5 shows a contour plot of the pressure on the fine 
grid. The regular pattern of the TS waves is clearly 
visible. The initial amplitude is very small and de- 
cays, because the disturbance is in the stable part of 
the boundary layer. More downstream the disturbance 
grows very regularly, saturates and decays. The decay 
is partly physical and at the end of the plate further 
increased by the coarsening of the grid to minimize re- 
flections from the outflow boundary. 

Figures 6 and 7 show the spatial growth rates obtained 
with the finite difference and finite volumes scheme on 
the fine grid and compare them with the PSE results of 
Bertolotti 2 . They are virtually identical till the distur- 
bances reach the buffer region, where the growth rate 
suddenly changes and the Tollmien-Schlichting wave 
rapidly decays. The CPU time used for both schemes 
was approximately equal for the implicit calculations, 
with the finite volume scheme 1.1 times more expensive 
than the finite difference scheme. The approximation of 
the flux integral in the finite volume method, equation 
(3.12), however, is crucial to minimize computing time 
for the finite volume scheme. These plots also show 
that the buffer region is quite effective in minimizing 
reflections from the outer wall. The simulation was 
continued in all cases till the leading wave front would 
have travelled at least twice the length of the domain. 
The sound waves, which travel faster, then would have 
a chance to reflect several times through the domain, 
with no apparent effect on the growth rates. The use 
of a buffer region is, however, not without pitfalls. One 
has to be very careful to create a smooth transition 
with the flow domain. 

The accuracy of the simulation also depends on the 
time step and time integration scheme used. The time 
integration method is a second order 3 point implicit 
multi-step scheme. Four Newton iterations were used 
to improve time accuracy of the implicit scheme. The 
residue decreased two orders of magnitude during the 
Newton iterations and was approximately 5 x 10" 7 at 
the end of the Newton iterations. Figures 8 and 9 show 


the growth rate of the normal and streamwise velocity 
component for simulations with different time steps us- 
ing the finite volume scheme, but all on the 436 x 80 
grid. The time step At is equivalent with a CFL num- 
ber 80. As can be seen from these plots the accuracy 
is bounded by the time integration scheme and not by 
the spatial resolution. Significant improvement should 
be obtained by using a higher order accurate time in- 
tegration method. 

The results of computations with the finite difference 
scheme on the coarse and fine grid are presented in Fig- 
ure 10. They show that the coarse grid simulation is 
underresolved, while the fine grid simulation compares 
well with the PSE results of Bertolotti 2 . It should be 
emphasized that it is very important to perform the 
computations on different grid levels to test accuracy, 
especially when there are no theoretical results avail- 
able. The coarse grid results do show that the flow 
becomes unstable but the growth rate is not correct 
and can only be checked by increasing the resolution. 

Finally, the effect of the location of suction and blow- 
ing was investigated. Figure 11. shows results of sim- 
ulations with the finite difference scheme with suction 
and blowing applied directly after the inflow boundary 
or just before the region where the flow becomes un- 
stable. It can be seen that the growth rates are not 
sensitive to the location of suction and blowing. 

To summarize the results discussed in this section it 
can be stated that both the finite volume and finite dif- 
ference schemes can be used for simulations of bound- 
ary layer instability. The results for the incompressible 
flow were slightly better for the finite difference scheme, 
but this can be attributed to the second order viscous 
contribution in the finite volume scheme. The biggest 
advantage of the finite difference scheme is that it is less 
sensitive to grid stretching, but this scheme is not eas- 
ily extended to flows with sonic points, because it will 
loose accuracy at these points, which is not the case 
for the finite volume scheme. It was also found that 
using a higher order accurate time integration method 
will be more efficient and extreme care has to be taken 
to guarantee numerical accuracy, preferably by using 
grids with different resolution. 
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Figure 1. Stream wise velocity U at 10 equally Figure 2. Normal velocity V ^Re x at 10 equally 

spaced stations, Re x = 50.000 - 320.000, = spaced stations, Re x - 50.000 - 320.000, - 

.08, versus similarity parameter 77 = ^ yf Re, com- *08, versus similarity parameter 77 = ^ >J Re, com- 
pared with compressible Blasius solution (dashed pared with compressible Blasius solution (dashed 

line). line). 



Figure 3. Spatial growth rate — a* of streamwise 
velocity disturbances compared with results from 
Fasel et al. 5 . Moo = -08 
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Figure 4. Spatial growth rate — a* of normal 
velocity disturbances compared with results from 
Fasel et al. 5 . M^ — .08 
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Figure 5, Pressure contours in flat plate boundary layer, M » = .5, initial amplitude suction and blowing at wall 
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Figure 6. Spatial growth rate — a* of stream- 
wise velocity disturbances using finite volume and 
finite difference schemes compared with PSE re- 
sults Bertolotti 2 , = .5 


Figure 7. Spatial growth rate — ai of normal 
velocity disturbances using finite volume and fi- 
nite difference schemes compared with PSE results 
Bertolotti 2 , Moo — *5 
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Figure 8. Spatial growth rate —a* of streamwi se 
velocity disturbances using finite volume scheme 
with different time steps, Moo = .5 
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Figure 10. Spatial growth rate — of stream wise 
velocity disturbances using finite difference scheme 
on coarse and fine grid, — .5 
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Figure 9. Spatial growth rate -a* of normal 
velocity disturbances using finite volume scheme 
with different time steps, Moo — -5 
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Figure 11. Spatial growth rate — a* of stream wise 
velocity disturbances using finite difference scheme 
for different locations wall suction-blowing, Moo = 
.5 
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